---

In this notebook, we will show how to create pcraster maps to represent soil and land hydrologic properties for a SPHY model input.

---

00. Load Packages

import os, shutil
import pcraster as pcr
import pcr_tool

import numpy as np
import cartopy

01. Basic Setup for a Bhutan watershed (Punatsangchhu, Wangchhu, Amochhu, and Jaldakah)

epsg_model = 5266 # Bhutan National Grid: a transverse Mercator projection
epsg_lonlat = 4326 # the common lat-lon coordinate reference system

path_data = os.path.join('..', 'data')
path_output = os.path.join('..', 'input') # the "input" folder will be directly used as an input path for the SPHY model run
os.makedirs(path_output, exist_ok = True)

# plotting options
flag_plot_aguila = True
flag_plot_cartopy = True

# Bhutan Administrative Boundaries (Level 1) to be added in a plot
path_shape = os.path.join(path_data, 'SHAPE', 'btn_admbnda_adm1_bnlc_20201026.shp')
shape_geometry = cartopy.io.shapereader.Reader(path_shape).geometries()
shape_feature = cartopy.feature.ShapelyFeature(shape_geometry, crs = cartopy.crs.PlateCarree(), facecolor = 'none', linewidth = 0.25)

02. Load Clonemap and DEM

First, we need to declare the class that we will be using throughout this notebook, with an EPSG to define a clonemap.

print(f'- Input EPSG:\n{epsg_model}')
map_clone = pcr_tool.clonemap(epsg = epsg_model)

- Input EPSG:
5266

Let's load the clonemap and DEM that we generated in the previous notebook.

# load clonemap
path_clonemap = os.path.join(path_output, 'clone.map')
print(f'- Input PCRaster Path for Clonemap:\n{path_clonemap}')
map_clone.load_clonemap(path = path_clonemap)

# load DEM
path_dem = os.path.join(path_output, 'dem.map')
print(f'- Input PCRaster Path for DEM:\n{path_dem}')
map_dem = map_clone.read_pcraster(path = path_dem)
map_clone.set_dem(map_dem)

print('- Outputs:')
print(map_dem)

# plots
if flag_plot_aguila:
    pcr.aguila(map_clone.clonemap)
    pcr.aguila(map_dem)

- Input PCRaster Path for Clonemap:
../input/clone.map
- Input PCRaster Path for DEM:
../input/dem.map
- Outputs:
<pcraster._pcraster.Field object at 0x16c73e0f0>

03. Set Soil-Related Hydrologic Properties

First, we need to define the depths for the rootzone and subsoil layers. Let's assume the uniform depths for both layers.

rootdepth = 300 # rootzone depth (mm)
subdepth = 1700 # subsoil depth (mm)

array_rootdepth = rootdepth * np.ones([map_clone.nrow, map_clone.ncol]) # create numpy.ndarray with the same numbers of clonemap rows and columns
array_subdepth = subdepth * np.ones([map_clone.nrow, map_clone.ncol]) # create numpy.ndarray with the same numbers of clonemap rows and columns
# conversion from numpy.ndarray to pcraster.Field
map_rootdepth = map_clone.from_numpy(
    array = array_rootdepth, # numpy.ndarray to convert into pcraster.Field
    datatype = pcr.Scalar, # pcraster data type
    missing_value = 0, # missing value to assign
    ) # for rootzone layer
map_subdepth = map_clone.from_numpy(
    array = array_subdepth, # numpy.ndarray to convert into pcraster.Field
    datatype = pcr.Scalar, # pcraster data type
    missing_value = 0, # missing value to assign
    ) # for subsoil layer

print('- Outputs:')
print(map_rootdepth)
print(map_subdepth)

# export into a pcraster file
map_clone.write_pcraster(
    path = os.path.join(path_output, 'root_depth.map'), # path to write in the pcraster format
    field = map_rootdepth, # pcraster.Field to write
    )
map_clone.write_pcraster(
    path = os.path.join(path_output, 'sub_depth.map'), # path to write in the pcraster format
    field = map_subdepth, # pcraster.Field to write
    )

# plots
if flag_plot_aguila:
    pcr.aguila(map_rootdepth)
    pcr.aguila(map_subdepth)
if flag_plot_cartopy:
    map_clone.plot_cartopy(
        field = map_rootdepth, # pcraster.Field to plot
        shapes = [shape_feature], # a list of cartopy.feature.ShapelyFeature to plot together
        figsize = (6, 6), # matplotlib figure size
        basemap = cartopy.io.img_tiles.GoogleTiles(style = 'satellite'), # cartopy basemap (here, Google Satellite Images)
        basemap_level = 10, # cartopy basemap zoom level
        extent = None, # cartopy map extent
        vmin = None, # matplotlib contour min
        vmax = None, # matplotlib contour max
        cmap = 'turbo', # matplotlib contour colormap
        alpha = 0.7, # matplotlib contour transparency
        title = 'root_depth.map', # matplotlib figure title
        clabel = 'rootzone depth', # matplotlib contour label
        savefig = None, # path to save the figure
        )
    map_clone.plot_cartopy(
        field = map_subdepth, # pcraster.Field to plot
        shapes = [shape_feature], # a list of cartopy.feature.ShapelyFeature to plot together
        figsize = (6, 6), # matplotlib figure size
        basemap = cartopy.io.img_tiles.GoogleTiles(style = 'satellite'), # cartopy basemap (here, Google Satellite Images)
        basemap_level = 10, # cartopy basemap zoom level
        extent = None, # cartopy map extent
        vmin = None, # matplotlib contour min
        vmax = None, # matplotlib contour max
        cmap = 'turbo', # matplotlib contour colormap
        alpha = 0.7, # matplotlib contour transparency
        title = 'sub_depth.map', # matplotlib figure title
        clabel = 'subsoil depth', # matplotlib contour label
        savefig = None, # path to save the figure
        )

- Outputs:
<pcraster._pcraster.Field object at 0x105eae030>
<pcraster._pcraster.Field object at 0x16c377330>
write_pcraster: written to ../input/root_depth.map.
write_pcraster: written to ../input/sub_depth.map.

We're first going to generate soil-related hydrologic properties by assigning the values on the lookup table to the soil classifications (e.g., Harmonized World Soil Database (HWSD) V2).

path_soil = os.path.join(path_data, 'SOIL', 'hwsd2_subset.tif')

print(f'- Input GEOTIFF Path for Soil Classifications:\n{path_soil}')
# pass "epsg = epsg_lonlat" as the DEM file has a different CRS (lat/lon)
map_soil = map_clone.read_raster(
    path = path_soil, # geotiff path to read
    epsg = epsg_lonlat, # input EPSG
    datatype = pcr.Nominal, # pcraster data type
    band = 1, # band number to read
    ) # automatically subset the input based on the clonemap extent

print('- Outputs:')
print(map_soil)

# export into a pcraster file
map_clone.write_pcraster(
    path = os.path.join(path_output, 'soil.map'), # path to write in the pcraster format
    field = map_soil, # pcraster.Field to write
    )

# plots
if flag_plot_aguila: pcr.aguila(map_soil)
if flag_plot_cartopy:
    map_clone.plot_cartopy(
        field = map_soil, # pcraster.Field to plot
        shapes = [shape_feature], # a list of cartopy.feature.ShapelyFeature to plot together
        figsize = (6, 6), # matplotlib figure size
        basemap = cartopy.io.img_tiles.GoogleTiles(style = 'satellite'), # cartopy basemap (here, Google Satellite Images)
        basemap_level = 10, # cartopy basemap zoom level
        extent = None, # cartopy map extent
        vmin = None, # matplotlib contour min
        vmax = None, # matplotlib contour max
        cmap = 'turbo', # matplotlib contour colormap
        alpha = 0.7, # matplotlib contour transparency
        title = 'soil.map', # matplotlib figure title
        clabel = 'soil (HWSD V2)', # matplotlib contour label
        savefig = None, # path to save the figure
        )

- Input GEOTIFF Path for Soil Classifications:
../data/SOIL/hwsd2_subset.tif
- Outputs:
<pcraster._pcraster.Field object at 0x168b258f0>
write_pcraster: written to ../input/soil.map.

Currently, the hydrologic properties on the lookup tables are based on Harmonized World Soil Database (HWSD) V1, which requires mapping from HWSD V2 to HWSD V1.

path_table_HWSD2to1 = os.path.join(path_data, 'SOIL', 'HWSD2to1.tbl')

print(f'- Input Lookup Table Path for HWSD Mapping:\n{path_table_HWSD2to1}')
map_soil_ref = pcr.lookupscalar(path_table_HWSD2to1, map_soil) # map pcraster.Field from HWSD V2 to V1

print('- Outputs:')
print(map_soil_ref)

# export into a pcraster file
map_clone.write_pcraster(
    path = os.path.join(path_output, 'soilref.map'), # path to write in the pcraster format
    field = map_soil_ref # pcraster.Field to write
    )

#plots
if flag_plot_aguila: pcr.aguila(map_soil_ref)
if flag_plot_cartopy: map_clone.plot_cartopy(
    field = map_soil_ref, # pcraster.Field to plot
    shapes = [shape_feature], # a list of cartopy.feature.ShapelyFeature to plot together
    figsize = (6, 6), # matplotlib figure size
    basemap = cartopy.io.img_tiles.GoogleTiles(style = 'satellite'), # cartopy basemap (here, Google Satellite Images)
    basemap_level = 10, # cartopy basemap zoom level
    extent = None, # cartopy map extent
    vmin = None, # matplotlib contour min
    vmax = None, # matplotlib contour max
    cmap = 'turbo', # matplotlib contour colormap
    alpha = 0.7, # matplotlib contour transparency
    title = 'soilref.map', # matplotlib figure title
    clabel = 'Soil (HWSD V1.2)', # matplotlib contour label
    savefig = None, # path to save the figure
    )

- Input Lookup Table Path for HWSD Mapping:
../data/SOIL/HWSD2to1.tbl
- Outputs:
<pcraster._pcraster.Field object at 0x173f6c1f0>
write_pcraster: written to ../input/soilref.map.

Now, we can use the Harmonized World Soil Database (HWSD) V1 soil classifications to assign hydrologic properties from the HWSD V1 lookup tables. Let's generate soil-related hydrologic properties for rootzone and subsoil layers, respectively, based on the lookup tables defined for Harmonized World Soil Database (HWSD) V1.

• field: field capacity [mm/mm]
• pwlt: permanent wilting point [mm/mm]
• satk: saturated hydraulic conductivity [mm/day]
• satwc: saturated water content [mm/mm]
• wlt: wilting point [mm/mm]

# field capacity in rootzone
path_table_rfield = os.path.join(path_data, 'SOIL', 'root_field.tbl') # lookup table for field capacity in rootzone [mm/mm]
print(f'- Input Lookup Table Path for Rootzone Field Capacity:\n{path_table_rfield}')
map_rfield = pcr.lookupscalar(path_table_rfield, map_soil_ref) # assign the field capacity in rootzone from the corresponding lookup table

print('- Outputs:')
print(map_rfield)

# export into a pcraster file
map_clone.write_pcraster(
    path = os.path.join(path_output, 'root_field.map'), # path to write in the pcraster format
    field = map_rfield, # pcraster.Field to write
    )

# plots
if flag_plot_aguila: pcr.aguila(map_rfield)
if flag_plot_cartopy:
    map_clone.plot_cartopy(
        field = map_rfield, # pcraster.Field to plot
        shapes = [shape_feature], # a list of cartopy.feature.ShapelyFeature to plot together
        figsize = (6, 6), # matplotlib figure size
        basemap = cartopy.io.img_tiles.GoogleTiles(style = 'satellite'), # cartopy basemap (here, Google Satellite Images)
        basemap_level = 10, # cartopy basemap zoom level
        extent = None, # cartopy map extent
        vmin = None, # matplotlib contour min
        vmax = None, # matplotlib contour max
        cmap = 'turbo', # matplotlib contour colormap
        alpha = 0.7, # matplotlib contour transparency
        title = 'root_field.map', # matplotlib figure title
        clabel = 'Field Capacity in Rootzone [mm/mm]', # matplotlib contour label
        savefig = None, # path to save the figure
        )

- Input Lookup Table Path for Rootzone Field Capacity:
../data/SOIL/root_field.tbl
- Outputs:
<pcraster._pcraster.Field object at 0x173f52ab0>
write_pcraster: written to ../input/root_field.map.

# permanent wilting point in rootzone
path_table_rpwlt = os.path.join(path_data, 'SOIL', 'root_pwlt.tbl') # lookup table for permanent wilting point in rootzone [mm/mm]
print(f'- Input Lookup Table Path for Rootzone Permanent Wilting Point:\n{path_table_rpwlt}')
map_rpwlt = pcr.lookupscalar(path_table_rpwlt, map_soil_ref) # assign the permanent wilting point in rootzone on the corresponding lookup table

print('- Outputs:')
print(map_rpwlt)

# export into a pcraster file
map_clone.write_pcraster(
    path = os.path.join(path_output, 'root_pwlt.map'), # path to write in the pcraster format
    field = map_rpwlt, # pcraster.Field to write
    )

# plots
if flag_plot_aguila: pcr.aguila(map_rpwlt)
if flag_plot_cartopy:
    map_clone.plot_cartopy(
        field = map_rpwlt, # pcraster.Field to plot
        shapes = [shape_feature], # a list of cartopy.feature.ShapelyFeature to plot together
        figsize = (6, 6), # matplotlib figure size
        basemap = cartopy.io.img_tiles.GoogleTiles(style = 'satellite'), # cartopy basemap (here, Google Satellite Images)
        basemap_level = 10, # cartopy basemap zoom level
        extent = None, # cartopy map extent
        vmin = None, # matplotlib contour min
        vmax = None, # matplotlib contour max
        cmap = 'turbo', # matplotlib contour colormap
        alpha = 0.7, # matplotlib contour transparency
        title = 'root_pwlt.map', # matplotlib figure title
        clabel = 'Permanent Wilting Point in Rootzone [mm/mm]', # matplotlib contour label
        savefig = None, # path to save the figure
        )

- Input Lookup Table Path for Rootzone Permanent Wilting Point:
../data/SOIL/root_pwlt.tbl
- Outputs:
<pcraster._pcraster.Field object at 0x173f1abb0>
write_pcraster: written to ../input/root_pwlt.map.

# saturated hydraulic conductivity in rootzone
path_table_rsatk = os.path.join(path_data, 'SOIL', 'root_satk.tbl') # lookup table for saturated hydraulic conductivity in rootzone [mm/day]
print(f'- Input Lookup Table Path for Rootzone Saturated Hydraulic Conductivity:\n{path_table_rsatk}')
map_rsatk = pcr.lookupscalar(path_table_rsatk, map_soil_ref) # assign the saturated hydraulic conductivity in rootzone on the corresponding lookup table

print('- Outputs:')
print(map_rsatk)

# export into a pcraster file
map_clone.write_pcraster(
    path = os.path.join(path_output, 'root_satk.map'), # path to write in the pcraster format
    field = map_rsatk, # pcraster.Field to write
    )

# plots
if flag_plot_aguila: pcr.aguila(map_rsatk)
if flag_plot_cartopy:
    map_clone.plot_cartopy(
        field = map_rsatk, # pcraster.Field to plot
        shapes = [shape_feature], # a list of cartopy.feature.ShapelyFeature to plot together
        figsize = (6, 6), # matplotlib figure size
        basemap = cartopy.io.img_tiles.GoogleTiles(style = 'satellite'), # cartopy basemap (here, Google Satellite Images)
        basemap_level = 10, # cartopy basemap zoom level
        extent = None, # cartopy map extent
        vmin = None, # matplotlib contour min
        vmax = None, # matplotlib contour max
        cmap = 'turbo', # matplotlib contour colormap
        alpha = 0.7, # matplotlib contour transparency
        title = 'root_satk.map', # matplotlib figure title
        clabel = 'Saturated Hydraulic Conductivity in Rootzone [mm/day]', # matplotlib contour label
        savefig = None, # path to save the figure
        )

- Input Lookup Table Path for Rootzone Saturated Hydraulic Conductivity:
../data/SOIL/root_satk.tbl
- Outputs:
<pcraster._pcraster.Field object at 0x105edbcf0>
write_pcraster: written to ../input/root_satk.map.

# saturated water content in rootzone
path_table_rsatwc = os.path.join(path_data, 'SOIL', 'root_satwc.tbl') # lookup table for saturated water content in rootzone [mm/mm]
print(f'- Input Lookup Table Path for Rootzone Saturated Water Content:\n{path_table_rsatwc}')
map_rsatwc = pcr.lookupscalar(path_table_rsatwc, map_soil_ref) # assign the saturated water content in rootzone on the corresponding lookup table

print('- Outputs:')
print(map_rsatwc)

# export into a pcraster file
map_clone.write_pcraster(
    path = os.path.join(path_output, 'root_satwc.map'), # path to write in the pcraster format
    field = map_rsatwc, # pcraster.Field to write
    )

# plots
if flag_plot_aguila: pcr.aguila(map_rsatwc)
if flag_plot_cartopy:
    map_clone.plot_cartopy(
        field = map_rsatwc, # pcraster.Field to plot
        shapes = [shape_feature], # a list of cartopy.feature.ShapelyFeature to plot together
        figsize = (6, 6), # matplotlib figure size
        basemap = cartopy.io.img_tiles.GoogleTiles(style = 'satellite'), # cartopy basemap (here, Google Satellite Images)
        basemap_level = 10, # cartopy basemap zoom level
        extent = None, # cartopy map extent
        vmin = None, # matplotlib contour min
        vmax = None, # matplotlib contour max
        cmap = 'turbo', # matplotlib contour colormap
        alpha = 0.7, # matplotlib contour transparency
        title = 'root_satwc.map', # matplotlib figure title
        clabel = 'Saturated Water Content in Rootzone [mm/mm]', # matplotlib contour label
        savefig = None, # path to save the figure
        )

- Input Lookup Table Path for Rootzone Saturated Water Content:
../data/SOIL/root_satwc.tbl
- Outputs:
<pcraster._pcraster.Field object at 0x1755ea630>
write_pcraster: written to ../input/root_satwc.map.

# wilting point in rootzone
path_table_rwlt = os.path.join(path_data, 'SOIL', 'root_wlt.tbl') # lookup table for wilting point in rootzone [mm/mm]
print(f'- Input Lookup Table Path for Rootzone Wilting Point:\n{path_table_rwlt}')
map_rwlt = pcr.lookupscalar(path_table_rwlt, map_soil_ref) # assign the wilting point in rootzone on the corresponding lookup table

print('- Outputs:')
print(map_rwlt)

# export into a pcraster file
map_clone.write_pcraster(
    path = os.path.join(path_output, 'root_wlt.map'), # path to write in the pcraster format
    field = map_rwlt, # pcraster.Field to write
    )

# plots
if flag_plot_aguila: pcr.aguila(map_rwlt)
if flag_plot_cartopy:
    map_clone.plot_cartopy(
        field = map_rwlt, # pcraster.Field to plot
        shapes = [shape_feature], # a list of cartopy.feature.ShapelyFeature to plot together
        figsize = (6, 6), # matplotlib figure size
        basemap = cartopy.io.img_tiles.GoogleTiles(style = 'satellite'), # cartopy basemap (here, Google Satellite Images)
        basemap_level = 10, # cartopy basemap zoom level
        extent = None, # cartopy map extent
        vmin = None, # matplotlib contour min
        vmax = None, # matplotlib contour max
        cmap = 'turbo', # matplotlib contour colormap
        alpha = 0.7, # matplotlib contour transparency
        title = 'root_wlt.map', # matplotlib figure title
        clabel = 'Wilting Point in Rootzone [mm/mm]', # matplotlib contour label
        savefig = None, # path to save the figure
        )

- Input Lookup Table Path for Rootzone Wilting Point:
../data/SOIL/root_wlt.tbl
- Outputs:
<pcraster._pcraster.Field object at 0x1757067b0>
write_pcraster: written to ../input/root_wlt.map.

# field capacity in subsoil
path_table_sfield = os.path.join(path_data, 'SOIL', 'sub_field.tbl') # lookup table for field capacity in subsoil [mm/mm]
print(f'- Input Lookup Table Path for Subsoil Field Capacity:\n{path_table_sfield}')
map_sfield = pcr.lookupscalar(path_table_sfield, map_soil_ref) # assign the field capacity in subsoil on the corresponding lookup table

print('- Outputs:')
print(map_sfield)

# export into a pcraster file
map_clone.write_pcraster(
    path = os.path.join(path_output, 'sub_field.map'), # path to write in the pcraster format
    field = map_sfield, # pcraster.Field to write
    )

# plots
if flag_plot_aguila: pcr.aguila(map_sfield)
if flag_plot_cartopy:
    map_clone.plot_cartopy(
        field = map_sfield, # pcraster.Field to plot
        shapes = [shape_feature], # a list of cartopy.feature.ShapelyFeature to plot together
        figsize = (6, 6), # matplotlib figure size
        basemap = cartopy.io.img_tiles.GoogleTiles(style = 'satellite'), # cartopy basemap (here, Google Satellite Images)
        basemap_level = 10, # cartopy basemap zoom level
        extent = None, # cartopy map extent
        vmin = None, # matplotlib contour min
        vmax = None, # matplotlib contour max
        cmap = 'turbo', # matplotlib contour colormap
        alpha = 0.7, # matplotlib contour transparency
        title = 'sub_field.map', # matplotlib figure title
        clabel = 'Field Capacity in Subsoil [mm/mm]', # matplotlib contour label
        savefig = None, # path to save the figure
        )

- Input Lookup Table Path for Subsoil Field Capacity:
../data/SOIL/sub_field.tbl
- Outputs:
<pcraster._pcraster.Field object at 0x17556e730>
write_pcraster: written to ../input/sub_field.map.

# saturated hydraulic conductivity in subsoil
path_table_ssatk = os.path.join(path_data, 'SOIL', 'sub_satk.tbl') # lookup table for saturated hydraulic conductivity in subsoil [mm/day]
print(f'- Input Lookup Table Path for Subsoil Saturated Hydraulic Conductivity:\n{path_table_ssatk}')
map_ssatk = pcr.lookupscalar(path_table_ssatk, map_soil_ref) # assign the saturated hydraulic conductivity in subsoil on the corresponding lookup table

print('- Outputs:')
print(map_ssatk)

# export into a pcraster file
map_clone.write_pcraster(
    path = os.path.join(path_output, 'sub_satk.map'), # path to write in the pcraster format
    field = map_ssatk, # pcraster.Field to write
    )

# plots
if flag_plot_aguila: pcr.aguila(map_ssatk)
if flag_plot_cartopy:
    map_clone.plot_cartopy(
        field = map_ssatk, # pcraster.Field to plot
        shapes = [shape_feature], # a list of cartopy.feature.ShapelyFeature to plot together
        figsize = (6, 6), # matplotlib figure size
        basemap = cartopy.io.img_tiles.GoogleTiles(style = 'satellite'), # cartopy basemap (here, Google Satellite Images)
        basemap_level = 10, # cartopy basemap zoom level
        extent = None, # cartopy map extent
        vmin = None, # matplotlib contour min
        vmax = None, # matplotlib contour max
        cmap = 'turbo', # matplotlib contour colormap
        alpha = 0.7, # matplotlib contour transparency
        title = 'sub_satk.map', # matplotlib figure title
        clabel = 'Saturated Hydraulic Conductivity in Subsoil [mm/day]', # matplotlib contour label
        savefig = None, # path to save the figure
        )

- Input Lookup Table Path for Subsoil Saturated Hydraulic Conductivity:
../data/SOIL/sub_satk.tbl
- Outputs:
<pcraster._pcraster.Field object at 0x175610c70>
write_pcraster: written to ../input/sub_satk.map.

# saturated water content in subsoil
path_table_ssatwc = os.path.join(path_data, 'SOIL', 'sub_satwc.tbl') # lookup table for saturated water content in subsoil [mm/mm]
print(f'- Input Lookup Table Path for Subsoil Saturated Water Content:\n{path_table_ssatwc}')
map_ssatwc = pcr.lookupscalar(path_table_ssatwc, map_soil_ref) # assign the saturated water content in subsoil on the corresponding lookup table

print('- Outputs:')
print(map_ssatwc)

# export into a pcraster file
map_clone.write_pcraster(
    path = os.path.join(path_output, 'sub_satwc.map'), # path to write in the pcraster format
    field = map_ssatwc, # pcraster.Field to write
    )

# plots
if flag_plot_aguila: pcr.aguila(map_ssatwc)
if flag_plot_cartopy:
    map_clone.plot_cartopy(
        field = map_ssatwc, # pcraster.Field to plot
        shapes = [shape_feature], # a list of cartopy.feature.ShapelyFeature to plot together
        figsize = (6, 6), # matplotlib figure size
        basemap = cartopy.io.img_tiles.GoogleTiles(style = 'satellite'), # cartopy basemap (here, Google Satellite Images)
        basemap_level = 10, # cartopy basemap zoom level
        extent = None, # cartopy map extent
        vmin = None, # matplotlib contour min
        vmax = None, # matplotlib contour max
        cmap = 'turbo', # matplotlib contour colormap
        alpha = 0.7, # matplotlib contour transparency
        title = 'sub_satwc.map', # matplotlib figure title
        clabel = 'Saturated Water Content in Subsoil [mm/mm]', # matplotlib contour label
        savefig = None, # path to save the figure
        )

- Input Lookup Table Path for Subsoil Saturated Water Content:
../data/SOIL/sub_satwc.tbl
- Outputs:
<pcraster._pcraster.Field object at 0x1756ac1f0>
write_pcraster: written to ../input/sub_satwc.map.

ALTERNATIVE: If the hydrologic properties are available geospatially from a set of rasters (e.g.,HiHydroSoil 250m), we can directly generate them by importing the raster files.

# field capacity in rootzone
path_rfield_HHS = os.path.join(path_data, 'SOIL', 'WCpF2_M_250m_TOPSOIL_subset.tif') # raster for field capacity in rootzone [mm/mm]
print(f'- Input Raster Path for Rootzone Field Capacity:\n{path_rfield_HHS}')
# pass "epsg = epsg_check" as the input file has a different CRS (lat/lon)
map_rfield_HHS = map_clone.read_raster(
    path = path_rfield_HHS, # geotiff path to read
    epsg = epsg_lonlat, # input EPSG
    datatype = pcr.Scalar, # pcraster data type
    band = 1, # band number to read
    )

print('- Outputs:')
print(map_rfield_HHS)

# export into a pcraster file
map_clone.write_pcraster(
    path = os.path.join(path_output, 'root_field_HHS.map'), # path to write in the pcraster format
    field = map_rfield_HHS, # pcraster.Field to write
    )

# plots
if flag_plot_aguila: pcr.aguila(map_rfield_HHS)
if flag_plot_cartopy:
    map_clone.plot_cartopy(
        field = map_rfield_HHS, # pcraster.Field to plot
        shapes = [shape_feature], # a list of cartopy.feature.ShapelyFeature to plot together
        figsize = (6, 6), # matplotlib figure size
        basemap = cartopy.io.img_tiles.GoogleTiles(style = 'satellite'), # cartopy basemap (here, Google Satellite Images)
        basemap_level = 10, # cartopy basemap zoom level
        extent = None, # cartopy map extent
        vmin = None, # matplotlib contour min
        vmax = None, # matplotlib contour max
        cmap = 'turbo', # matplotlib contour colormap
        alpha = 0.7, # matplotlib contour transparency
        title = 'root_field_HHS.map', # matplotlib figure title
        clabel = 'Field Capacity in Rootzone [mm/mm]', # matplotlib contour label
        savefig = None, # path to save the figure
        )

- Input Raster Path for Rootzone Field Capacity:
../data/SOIL/WCpF2_M_250m_TOPSOIL_subset.tif
- Outputs:
<pcraster._pcraster.Field object at 0x176905b30>
write_pcraster: written to ../input/root_field_HHS.map.
<urlopen error [Errno 8] nodename nor servname provided, or not known>

# permanent wilting point in rootzone
path_rpwlt_HHS = os.path.join(path_data, 'SOIL', 'WCpF4.2_M_250m_TOPSOIL_subset.tif') # raster for permanent wilting point in rootzone [mm/mm]
print(f'- Input Raster Path for Rootzone Permanent Wilting Point:\n{path_rpwlt_HHS}')
# pass "epsg = epsg_check" as the input file has a different CRS (lat/lon)
map_rpwlt_HHS = map_clone.read_raster(
    path = path_rpwlt_HHS, # geotiff path to read
    epsg = epsg_lonlat, # input EPSG
    datatype = pcr.Scalar, # pcraster data type
    band = 1, # band number to read
    )

print('- Outputs:')
print(map_rpwlt_HHS)

# export into a pcraster file
map_clone.write_pcraster(
    path = os.path.join(path_output, 'root_pwlt_HHS.map'), # path to write in the pcraster format
    field = map_rpwlt_HHS, # pcraster.Field to write
    )

# plots
if flag_plot_aguila: pcr.aguila(map_rpwlt_HHS)
if flag_plot_cartopy:
    map_clone.plot_cartopy(
        field = map_rpwlt_HHS, # pcraster.Field to plot
        shapes = [shape_feature], # a list of cartopy.feature.ShapelyFeature to plot together
        figsize = (6, 6), # matplotlib figure size
        basemap = cartopy.io.img_tiles.GoogleTiles(style = 'satellite'), # cartopy basemap (here, Google Satellite Images)
        basemap_level = 10, # cartopy basemap zoom level
        extent = None, # cartopy map extent
        vmin = None, # matplotlib contour min
        vmax = None, # matplotlib contour max
        cmap = 'turbo', # matplotlib contour colormap
        alpha = 0.7, # matplotlib contour transparency
        title = 'root_pwlt_HHS.map', # matplotlib figure title
        clabel = 'Permanent Wilting Point in Rootzone [mm/mm]', # matplotlib contour label
        savefig = None, # path to save the figure
        )

- Input Raster Path for Rootzone Permanent Wilting Point:
../data/SOIL/WCpF4.2_M_250m_TOPSOIL_subset.tif
- Outputs:
<pcraster._pcraster.Field object at 0x16c704fb0>
write_pcraster: written to ../input/root_pwlt_HHS.map.

# saturated hydraulic conductivity in rootzone
path_rsatk_HHS = os.path.join(path_data, 'SOIL', 'Ksat_M_250m_TOPSOIL_subset.tif') # raster for saturated hydraulic conductivity in rootzone [mm/day]
print(f'- Input Raster Path for Rootzone Saturated Hydraulic Conductivity:\n{path_rsatk_HHS}')
# pass "epsg = epsg_check" as the input file has a different CRS (lat/lon)
map_rsatk_HHS = map_clone.read_raster(
    path = path_rsatk_HHS, # geotiff path to read
    epsg = epsg_lonlat, # input EPSG
    datatype = pcr.Scalar, # pcraster data type
    band = 1, # band number to read
    )

print('- Outputs:')
print(map_rsatk_HHS)

# export into a pcraster file
map_clone.write_pcraster(
    path = os.path.join(path_output, 'root_satk_HHS.map'), # path to write in the pcraster format
    field = map_rsatk_HHS, # pcraster.Field to write
    )

# plots
if flag_plot_aguila: pcr.aguila(map_rsatk_HHS)
if flag_plot_cartopy:
    map_clone.plot_cartopy(
        field = map_rsatk_HHS, # pcraster.Field to plot
        shapes = [shape_feature], # a list of cartopy.feature.ShapelyFeature to plot together
        figsize = (6, 6), # matplotlib figure size
        basemap = cartopy.io.img_tiles.GoogleTiles(style = 'satellite'), # cartopy basemap (here, Google Satellite Images)
        basemap_level = 10, # cartopy basemap zoom level
        extent = None, # cartopy map extent
        vmin = None, # matplotlib contour min
        vmax = None, # matplotlib contour max
        cmap = 'turbo', # matplotlib contour colormap
        alpha = 0.7, # matplotlib contour transparency
        title = 'root_satk_HHS.map', # matplotlib figure title
        clabel = 'Saturated Hydraulic Conductivity in Rootzone [mm/day]', # matplotlib contour label
        savefig = None, # path to save the figure
        )

- Input Raster Path for Rootzone Saturated Hydraulic Conductivity:
../data/SOIL/Ksat_M_250m_TOPSOIL_subset.tif
- Outputs:
<pcraster._pcraster.Field object at 0x16c7326b0>
write_pcraster: written to ../input/root_satk_HHS.map.

# saturated water content in rootzone
path_rsatwc_HHS = os.path.join(path_data, 'SOIL', 'WCsat_M_250m_TOPSOIL_subset.tif') # raster for saturated water content in rootzone [mm/mm]
print(f'- Input Raster Path for Rootzone Saturated Water Content:\n{path_rsatwc_HHS}')
# pass "epsg = epsg_check" as the input file has a different CRS (lat/lon)
map_rsatwc_HHS = map_clone.read_raster(
    path = path_rsatwc_HHS, # geotiff path to read
    epsg = epsg_lonlat, # input EPSG
    datatype = pcr.Scalar, # pcraster data type
    band = 1, # band number to read
    )

print('- Outputs:')
print(map_rsatwc_HHS)

# export into a pcraster file
map_clone.write_pcraster(
    path = os.path.join(path_output, 'root_satwc_HHS.map'), # path to write in the pcraster format
    field = map_rsatwc_HHS, # pcraster.Field to write
    )

# plots
if flag_plot_aguila: pcr.aguila(map_rsatwc_HHS)
if flag_plot_cartopy:
    map_clone.plot_cartopy(
        field = map_rsatwc_HHS, # pcraster.Field to plot
        shapes = [shape_feature], # a list of cartopy.feature.ShapelyFeature to plot together
        figsize = (6, 6), # matplotlib figure size
        basemap = cartopy.io.img_tiles.GoogleTiles(style = 'satellite'), # cartopy basemap (here, Google Satellite Images)
        basemap_level = 10, # cartopy basemap zoom level
        extent = None, # cartopy map extent
        vmin = None, # matplotlib contour min
        vmax = None, # matplotlib contour max
        cmap = 'turbo', # matplotlib contour colormap
        alpha = 0.7, # matplotlib contour transparency
        title = 'root_satwc_HHS.map', # matplotlib figure title
        clabel = 'Saturated Water Content in Rootzone [mm/mm]', # matplotlib contour label
        savefig = None, # path to save the figure
        )

- Input Raster Path for Rootzone Saturated Water Content:
../data/SOIL/WCsat_M_250m_TOPSOIL_subset.tif
- Outputs:
<pcraster._pcraster.Field object at 0x175a8d430>
write_pcraster: written to ../input/root_satwc_HHS.map.

# wilting point in rootzone
path_rwlt_HHS = os.path.join(path_data, 'SOIL', 'WCpF3_M_250m_TOPSOIL_subset.tif') # raster for wilting point in rootzone [mm/mm]
print(f'- Input Raster Path for Rootzone Wilting Point:\n{path_rwlt_HHS}')
# pass "epsg = epsg_check" as the input file has a different CRS (lat/lon)
map_rwlt_HHS = map_clone.read_raster(
    path = path_rwlt_HHS, # geotiff path to read
    epsg = epsg_lonlat, # input EPSG
    datatype = pcr.Scalar, # pcraster data type
    band = 1, # band number to read
    )

print('- Outputs:')
print(map_rwlt_HHS)

# export into a pcraster file
map_clone.write_pcraster(
    path = os.path.join(path_output, 'root_wlt_HHS.map'), # path to write in the pcraster format
    field = map_rwlt_HHS, # pcraster.Field to write
    )

# plots
if flag_plot_aguila: pcr.aguila(map_rwlt_HHS)
if flag_plot_cartopy:
    map_clone.plot_cartopy(
        field = map_rwlt_HHS, # pcraster.Field to plot
        shapes = [shape_feature], # a list of cartopy.feature.ShapelyFeature to plot together
        figsize = (6, 6), # matplotlib figure size
        basemap = cartopy.io.img_tiles.GoogleTiles(style = 'satellite'), # cartopy basemap (here, Google Satellite Images)
        basemap_level = 10, # cartopy basemap zoom level
        extent = None, # cartopy map extent
        vmin = None, # matplotlib contour min
        vmax = None, # matplotlib contour max
        cmap = 'turbo', # matplotlib contour colormap
        alpha = 0.7, # matplotlib contour transparency
        title = 'root_wlt_HHS.map', # matplotlib figure title
        clabel = 'Wilting Point in Rootzone [mm/mm]', # matplotlib contour label
        savefig = None, # path to save the figure
        )

- Input Raster Path for Rootzone Wilting Point:
../data/SOIL/WCpF3_M_250m_TOPSOIL_subset.tif
- Outputs:
<pcraster._pcraster.Field object at 0x17567ef30>
write_pcraster: written to ../input/root_wlt_HHS.map.
<urlopen error [Errno 8] nodename nor servname provided, or not known>

# field capacity in subsoil
path_sfield_HHS = os.path.join(path_data, 'SOIL', 'WCpF2_M_250m_SUBSOIL_subset.tif') # raster for field capacity in subsoil [mm/mm]
print(f'- Input Raster Path for Subsoil Field Capacity:\n{path_sfield_HHS}')
# pass "epsg = epsg_check" as the input file has a different CRS (lat/lon)
map_sfield_HHS = map_clone.read_raster(
    path = path_sfield_HHS, # geotiff path to read
    epsg = epsg_lonlat, # input EPSG
    datatype = pcr.Scalar, # pcraster data type
    band = 1, # band number to read
    )

print('- Outputs:')
print(map_sfield_HHS)

# export into a pcraster file
map_clone.write_pcraster(
    path = os.path.join(path_output, 'sub_field_HHS.map'), # path to write in the pcraster format
    field = map_sfield_HHS, # pcraster.Field to write
    )

# plots
if flag_plot_aguila: pcr.aguila(map_sfield_HHS)
if flag_plot_cartopy:
    map_clone.plot_cartopy(
        field = map_sfield_HHS, # pcraster.Field to plot
        shapes = [shape_feature], # a list of cartopy.feature.ShapelyFeature to plot together
        figsize = (6, 6), # matplotlib figure size
        basemap = cartopy.io.img_tiles.GoogleTiles(style = 'satellite'), # cartopy basemap (here, Google Satellite Images)
        basemap_level = 10, # cartopy basemap zoom level
        extent = None, # cartopy map extent
        vmin = None, # matplotlib contour min
        vmax = None, # matplotlib contour max
        cmap = 'turbo', # matplotlib contour colormap
        alpha = 0.7, # matplotlib contour transparency
        title = 'sub_field_HHS.map', # matplotlib figure title
        clabel = 'Field Capacity in Subsoil [mm/mm]', # matplotlib contour label
        savefig = None, # path to save the figure
        )

- Input Raster Path for Subsoil Field Capacity:
../data/SOIL/WCpF2_M_250m_SUBSOIL_subset.tif
- Outputs:
<pcraster._pcraster.Field object at 0x17753ca30>
write_pcraster: written to ../input/sub_field_HHS.map.

# saturated hydraulic conductivity in subsoil
path_ssatk_HHS = os.path.join(path_data, 'SOIL', 'Ksat_M_250m_SUBSOIL_subset.tif') # raster for saturated hydraulic conductivity in subsoil [mm/day]
print(f'- Input Raster Path for Subsoil Saturated Hydraulic Conductivity:\n{path_ssatk_HHS}')
# pass "epsg = epsg_check" as the input file has a different CRS (lat/lon)
map_ssatk_HHS = map_clone.read_raster(
    path = path_ssatk_HHS, # geotiff path to read
    epsg = epsg_lonlat, # input EPSG
    datatype = pcr.Scalar, # pcraster data type
    band = 1, # band number to read
    )

print('- Outputs:')
print(map_ssatk_HHS)

# export into a pcraster file
map_clone.write_pcraster(
    path = os.path.join(path_output, 'sub_satk_HHS.map'), # path to write in the pcraster format
    field = map_ssatk_HHS, # pcraster.Field to write
    )

# plots
if flag_plot_aguila: pcr.aguila(map_ssatk_HHS)
if flag_plot_cartopy:
    map_clone.plot_cartopy(
        field = map_ssatk_HHS, # pcraster.Field to plot
        shapes = [shape_feature], # a list of cartopy.feature.ShapelyFeature to plot together
        figsize = (6, 6), # matplotlib figure size
        basemap = cartopy.io.img_tiles.GoogleTiles(style = 'satellite'), # cartopy basemap (here, Google Satellite Images)
        basemap_level = 10, # cartopy basemap zoom level
        extent = None, # cartopy map extent
        vmin = None, # matplotlib contour min
        vmax = None, # matplotlib contour max
        cmap = 'turbo', # matplotlib contour colormap
        alpha = 0.7, # matplotlib contour transparency
        title = 'sub_satk_HHS.map', # matplotlib figure title
        clabel = 'Saturated Hydraulic Conductivity in Subsoil [mm/day]', # matplotlib contour label
        savefig = None, # path to save the figure
        )

- Input Raster Path for Subsoil Saturated Hydraulic Conductivity:
../data/SOIL/Ksat_M_250m_SUBSOIL_subset.tif
- Outputs:
<pcraster._pcraster.Field object at 0x177858fb0>
write_pcraster: written to ../input/sub_satk_HHS.map.

# saturated water content in subsoil
path_ssatwc_HHS = os.path.join(path_data, 'SOIL', 'WCsat_M_250m_SUBSOIL_subset.tif') # raster for saturated water content in subsoil [mm/mm]
print(f'- Input Raster Path for Subsoil Saturated Water Content:\n{path_ssatwc_HHS}')
# pass "epsg = epsg_check" as the input file has a different CRS (lat/lon)
map_ssatwc_HHS = map_clone.read_raster(
    path = path_ssatwc_HHS, # geotiff path to read
    epsg = epsg_lonlat, # input EPSG
    datatype = pcr.Scalar, # pcraster data type
    band = 1, # band number to read
    )

print('- Outputs:')
print(map_ssatwc_HHS)

# export into a pcraster file
map_clone.write_pcraster(
    path = os.path.join(path_output, 'sub_satwc_HHS.map'), # path to write in the pcraster format
    field = map_ssatwc_HHS, # pcraster.Field to write
    )

# plots
if flag_plot_aguila: pcr.aguila(map_ssatwc_HHS)
if flag_plot_cartopy:
    map_clone.plot_cartopy(
        field = map_ssatwc_HHS, # pcraster.Field to plot
        shapes = [shape_feature], # a list of cartopy.feature.ShapelyFeature to plot together
        figsize = (6, 6), # matplotlib figure size
        basemap = cartopy.io.img_tiles.GoogleTiles(style = 'satellite'), # cartopy basemap (here, Google Satellite Images)
        basemap_level = 10, # cartopy basemap zoom level
        extent = None, # cartopy map extent
        vmin = None, # matplotlib contour min
        vmax = None, # matplotlib contour max
        cmap = 'turbo', # matplotlib contour colormap
        alpha = 0.7, # matplotlib contour transparency
        title = 'sub_satwc_HHS.map', # matplotlib figure title
        clabel = 'Saturated Water Content in Subsoil [mm/mm]', # matplotlib contour label
        savefig = None, # path to save the figure
        )

- Input Raster Path for Subsoil Saturated Water Content:
../data/SOIL/WCsat_M_250m_SUBSOIL_subset.tif
- Outputs:
<pcraster._pcraster.Field object at 0x1779529b0>
write_pcraster: written to ../input/sub_satwc_HHS.map.

04. Set Land-Related Hydrologic Properties

Next, we're going to generate a land use/cover pcraster map (e.g., NASA MODIS Land Cover), which will be combined with a lookup table for the crop coefficients by the SPHY model.

path_luse = os.path.join(path_data, 'LAND', 'MCD12Q1.A2001001.tif')
print(f'- Input Raster Path for Land Use/Cover:\n{path_luse}')
# pass "epsg = epsg_check" as the input file has a different CRS (lat/lon)
map_luse = map_clone.read_raster(
    path = path_luse, # geotiff path to read
    epsg = epsg_lonlat, # input EPSG
    datatype = pcr.Nominal, # pcraster data type
    band = 1, # band number to read
    )

print('- Outputs:')
print(map_luse)

# export into a pcraster file
map_clone.write_pcraster(
    path = os.path.join(path_output, 'luse.map'), # path to write in the pcraster format
    field = map_luse, # pcraster.Field to write
    )

# plots
if flag_plot_aguila: pcr.aguila(map_luse)
if flag_plot_cartopy:
    map_clone.plot_cartopy(
        field = map_luse, # pcraster.Field to plot
        shapes = [shape_feature], # a list of cartopy.feature.ShapelyFeature to plot together
        figsize = (6, 6), # matplotlib figure size
        basemap = cartopy.io.img_tiles.GoogleTiles(style = 'satellite'), # cartopy basemap (here, Google Satellite Images)
        basemap_level = 10, # cartopy basemap zoom level
        extent = None, # cartopy map extent
        vmin = None, # matplotlib contour min
        vmax = None, # matplotlib contour max
        cmap = 'turbo', # matplotlib contour colormap
        alpha = 0.7, # matplotlib contour transparency
        title = 'luse.map', # matplotlib figure title
        clabel = 'Land Use', # matplotlib contour label
        savefig = None, # path to save the figure
        )

- Input Raster Path for Land Use/Cover:
../data/LAND/MCD12Q1.A2001001.tif
- Outputs:
<pcraster._pcraster.Field object at 0x17759c3b0>
write_pcraster: written to ../input/luse.map.

Then, we need to define the crop coefficients in kc.tbl with respect to the land use/cover classification. You will see that the crop coefficients are arbitrarily assigned in kc.tbl, which requires furture updates through estimations or calibrations.

with open(file = os.path.join(path_data, 'kc.tbl'), mode = 'r') as f:
    for l in f:
        print(l)

0 0.0

1 1.0

2 1.0

3 1.0

4 1.0

5 1.0

6 1.0

7 1.0

8 1.0

9 1.0

10 1.0

11 1.0

12 1.0

13 1.0

14 1.0

15 1.0

16 1.0

254 0.0

255 0.0

Anyway, please do not forget to copy kc.tbl into the SPHY input folder.

shutil.copyfile(src = os.path.join(path_data, 'kc.tbl'), dst = os.path.join(path_output, 'kc.tbl'))

'../input/kc.tbl'

Lastly, the SPHY model requires a latitude pcraster map to calculate extraterrestrial radiation.

map_latitude = map_clone.create_latitude()

print('- Outputs:')
print(map_latitude)

# export into a pcraster file
map_clone.write_pcraster(
    path = os.path.join(path_output, 'latitude.map'), # path to write in the pcraster format
    field = map_latitude, # pcraster.Field to write
    )

# plots
if flag_plot_aguila: pcr.aguila(map_latitude)
if flag_plot_cartopy:
    map_clone.plot_cartopy(
        field = map_latitude, # pcraster.Field to plot
        shapes = [shape_feature], # a list of cartopy.feature.ShapelyFeature to plot together
        figsize = (6, 6), # matplotlib figure size
        basemap = cartopy.io.img_tiles.GoogleTiles(style = 'satellite'), # cartopy basemap (here, Google Satellite Images)
        basemap_level = 10, # cartopy basemap zoom level
        extent = None, # cartopy map extent
        vmin = None, # matplotlib contour min
        vmax = None, # matplotlib contour max
        cmap = 'turbo', # matplotlib contour colormap
        alpha = 0.7, # matplotlib contour transparency
        title = 'latitude.map', # matplotlib figure title
        clabel = 'Latitude [degN]', # matplotlib contour label
        savefig = None, # path to save the figure
        )

- Outputs:
<pcraster._pcraster.Field object at 0x177994230>
write_pcraster: written to ../input/latitude.map.